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Oscillations in a Simple Equatorial Climate System 

Ka-Ming Lau 

Department of Meteorology , Naval Postgraduate School, Monterey , CA 93940 
(Manuscript received 4 September 1979, in final form 14 October 1980) 

ABSTRACT 

Coupling between large-scale atmospheric and oceanic equatorial Kelvin waves is shown to be relevant 
in the climatic time scale related to equatorial oceaa'atmosphere processes. The present analyses show that 
the inclusion of air-sea coupling into the linearized shallow-water equations can result in two types of 
dispersion relations for the Kelvin waves. The first type (Mode I) has fast phase speed and is mostly 
manifest in the atmospheric response. This mode is relatively unaffected by air-sea coupling. The second 
type (Mode II) has slow phase speed and is the predominant mode in the time-scale of the ocean. A resonant 
stationary wave is shown to exist as an intrinsic mode in the coupled system, the length scale of which is 
determined by the strength of the coupling and the magnitude of the damping. The positive feedback 
mechanism shown to exist between these coupled Kelvin waves in the Mode II regime is suggested to 
play an important role in relation to observed low-latitude teleconnections. 

Results of the numerical experiments using the coupled model show that an El Nino-type oscillation 
can occur in a baroclinic ocean-atmosphere system as a result of a prolonged period of strengthening and 
subsequent weakening in the barotropic component of the wind. The weakening and the eastward shift of 
the rising branch of the Walker cell, identified as the atmospheric component of a coupled quasi-stationary 
Kelvin wave, during El Nino provide a positive feedback favoring warm water formation in the eastern 
Pacific and contribute to the large amplitude of the oscillation. 



1. Introduction 

In the study of the interannual fluctuations in the 
ocean- atmosphere system, the El Nino and associated 
events have attracted growing interests and concern 
among meteorologists and oceanographers largely 
because of the scientific importance as well as the 
far-reaching economic ramifications that have been 
attached to the phenomenon. Accompanying the El 
Nino anomalous SST changes are numerous ob- 
served changes such as the weakening of the trade 
winds, intensification of the equatorial countercur- 
rents, rapid rising sea level in the eastern Pacific 
and falling sea level to the west (Ramage, 1975; 
Wyrtki, 1975, 1977). The interrelationships between 
these various parameters and the causes of their 
changes are widely debated (Barnett, 1977). At pres- 
ent the leading contender of the various theories 
regarding the origin of El Nino is perhaps the one 
due to Wyrtki (1975). He maintained that the sudden 
onset of El Nino is only a signal of the dynamical 
response of the entire equatorial Pacific Ocean to 
atmospheric forcings. Because of the observation 
that a prolonged period of increased surface east- 
erlies always occurs in the central and western Pacific 
prior to the onset of each major El Nino, he proposed 
that these increased easterlies increase the strength 
of the westward flowing South Equatorial Current 
which tends to pile up water in the warmer western 



Pacific and increases the east-west sea level slope. 
Following the relaxation of the surface easterlies, 
the pressure gradient from the enhanced sea slope 
will drive back water toward the east leading to a 
rise in sea level and deepening of the warm surface 
layer along the eastern border of the Pacific. 

In equatorial oceans, Godfrey (1975), Hurlbert 
et al. (1976) and McCreary (1976) showed that the 
eastward movement of the water mass as a result of a 
relaxed surface easterly wind stress is likely to occur 
in the form of an eastward propagating equatorial 
Kelvin wave originating from the western boundary. 
The phase speed of this wave is ~2 m s" 1 which 
means a time of approximately two months for the 
wave to travel across the span of the equatorial Pacific 
Ocean. Such a time scale is consistent with the 
time lag between the relaxation of wind in the west 
and first occurrence of warm water in the east ob- 
served during El Nino. Study of the baroclinic re- 
sponse of the ocean by Cane and Sarachik (1976) 
also confirmed that the internal Kelvin mode plays a 
decisive role in determining the initial response of a 
disturbance originating from the western part of a 
meridionally bounded ocean as a result of a sym- 
metric weakening of zonal winds. They showed that 
as a Kelvin wave moves from west to east along the 
equator, water mass is continually built up ahead of 
the wave and an equal amount depleted behind it. 
At the coastal boundaries the reflection of the waves 
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produces, among others, equatorial-trapped Rossby 
waves whose phase speeds are about three times 
smaller than that of the Kelvin wave. In the absence 
of other external forcings, a relaxation in wind stress 
will result in an oversimplified manner, in a seiche- 
like oscillation across the ocean basin made up of 
fast eastward traveling Kelvin waves and slower 
westward propagating Rossby waves. 

On the other hand, in the equatorial troposphere 
Webster (1973) showed that the Walker circulation 
may be viewed as a stationary wave set up by reso- 
nant excitation forced by the semi-permanent heat 
sources in the tropics, and this east-west overturn- 
ing favors the largest zonal scale. Chang (1977) 
further demonstrated the importance of viscous 
damping to produce the correct phase speed and 
vertical scales in the observed low-frequency Kelvin- 
like oscillations of the forced equatorial waves. The 
variabilities in the amplitude and positions of the 
Walker circulation are then manifest in the response 
of stationary waves to slow variations in heat sources 
and sinks controlled, in part, by the change in sur- 
face conditions over the ocean. 

There are many El Nino-like events in the ocean- 
atmosphere system which exhibit dynamic similari- 
ties to the above but with different intensities and 
time scales (Wyrtki et al . , 1976). In particular, the 
time scale of the theoretically predicted equatorial 
oceanic waves is of the order of months, while the 
observed time scale of the fluctuations is of the order 
of years. It is reasonable to suggest that the coupling 
of the ocean-atmosphere system may introduce addi- 
tional important time scales. In view of the largely 
independent observations and hypotheses described 
above, there is a great need for studies to bring 
together these different hypotheses so that a coherent 
picture may hopefully emerge. This paper represents 
a first step toward this goal by using a simple coupled 
two-component (ocean/atmosphere) model which 
includes most of the elements of the physical proc- 
esses discussed above including the interaction be- 
tween the two components. The study is aimed at 
the large-scale spatial and temporal response of the 
coupled system to free and forced conditions, with 
emphasis on seeking closed cycles, if any exist, un- 
der certain given causal sequences. 

In Section 2 the model equations are developed 
and in Sections 3 and 4 some properties of the analyti- 
cal solutions, in particular those related to the mech- 
anisms of the coupling are discussed. Sections 5 and 
6 are devoted to the results of several numerical 
experiments designed to further elucidate, in terms 
of the results of the previous sections, the role of 
the air-sea coupling, in particular those related to 
El Nino, Including the effect of diabatic heating. 
This study parallels a recent paper by Lau (1979) in 
which the changes in Hadley-type circulations in 
relation to local changes of wind stress and oceanic 



upwelling were studied using a domain-averaged 
air-sea interactive model. 



2. The model 

For the purpose of elucidating the basic mech- 
anism involved in the coupling of large scale oceanic 
and atmospheric Kelvin waves, we use the following 
simple model. 



a. The atmosphere 

Since the equatorial Kelvin waves are character- 
ized by a vanishing meridional velocity component, 
we may write the linearized equations governing the 
Kelvin mode for a layer of incompressible fluid of 
homogeneous density with a free surface under hy- 
drostatic balance, as 



OUg 

"aT 



d<t> ^ 

+ — + DgUg 

dx 



Fa , 



d(j) du„ 

+ + Da<t> = Q a , 

dt dx 



( 1 ) 

( 2 ) 



where u a and <fr are the perturbation zonal velocity 
and geopotential; H a the equivalent depth, D a the 
linear damping coefficients and F„, Q a the coupling 
and forcing functions to be specified. Although (1)- 
(2) represent the so-called barotropic equations, they 
can be interpreted as the baroclinic component of a 
two-layer system in which u a and </> are defined by 

tta ~ W a i* 0 = ^2 $1> (2) 

where subscripts 2 and 1 refer to the upper and lower 
layer, respectively. The phase speed of the internal 
gravity wave in such case is given by c a = ( gH a ) m 
= [S(Ap) 2 /2] 1/2 , depending on the static stability S 
and thickness A p between the two levels. Therefore, 
(1) and (2) are also valid for internal mode of motions 
provided the velocity is taken as a wind shear and the 
geopotential replaced by the thickness or the mean 
temperature of the layers (Matsuno, 1966). In fact, 
(l)-(2) can have a more general interpretation ap- 
plicable to continuous, stratified fluid. In such a case 
the equivalent depth H a would have arisen as the 
separation constant or eigenvalue of the vertical 
structure equation in equatorial /3-plane theory. 
However, because the emphasis of this paper is on 
the ocean-atmosphere coupling mechanism, the study 
of the vertical structure is not pursued. 



b. The ocean 

The governing equations for an ocean compatible 
with our simplistic atmosphere model are the linear 
shallow water equations (e.g., Lighthill, 1969). The 
one-dimensional form of these equations linearized 
about a state of rest and uniform stratification are 
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Fig. 1 . Schematic diagram showing the structure of the coupled ocean-atmosphere 
system. See text for definitions of symbols. " 



given by 



dh 

+ g + D 0 u 0 = F 09 

dx 


(4) 


+ H 0 ^2 + D 0 h = Q 0t 
dx 


(5) 



where u 0 denotes the zonal current, h the deviation 
of the thermocline depth from its mean value, g the 
acceleration of gravity, D 0 the linear damping 
coefficients. The forcing functions F 0 and Q 0 which 
appear on the right-hand sides of (4)-(5) are 
the projections of the atmospheric forcing into 
the baroclinic mode of the ocean with equivalent 
depth H 0 which is of the order of 0.5 m or less in the 
equatorial regions (Moore and Philander, 1977). Eqs. 
(4) and (5) describe the time dependence of the Kelvin 
wave mode with a nondispersive phase speed given 
by c 0 = (gH 0 ) m . In a two-layer formulation the vari- 
ables u and h can be identified with a baroclinic 
zonal current and a depth-averaged temperature, re- 
spectively (Godfrey, 1975), i.e., 

«o = wot “ "o,; h = a(D x 0 i + D 2 0 2 )> (6) 

where subscripts 1, 2 refer to the top and bottom 
layers, respectively, d is the departure of the layer 
temperature from some mean value, a the coefficient 
of thermal expansion of water and D the depth of 
the oceanic layer. The dual nature of the thermocline 
depth /i as a measure of pressure (gh) and depth- 
averaged temperature [by (6)] allows a simple physi- 
cal interpretation of the response of the model in 
relation to the real ocean. In this paper we shall 
refer to h as the thermocline depth, sea level or 
SST where appropriate. The schematic structure of 
the coupled model is shown in Fig. 1. 

3. Solution of the coupled equations 

In order to solve the coupled shallow- water equa- 
tions (1), (2), (4) and (5), the first task is to specify 



the forcing functions and coupling terms. We assume 
the following simple forms, i.e.. 



F a = 0, 



Qa = a'^ a h + Q c f 



F 0 = — Q 0 = 0 

PoDi 



(7) 



where | V* | is a typical wind speed at the ocean sur- 
face; p a , p 0 the density of surface air and water, 
respectively; C D the drag coefficient, u s the zonal 
wind at the surface. The symbol represents the 
projection of the oceanic Kelvin mode onto the at- 
mosphere Kelvin mode and d>% vice versa. Here, 
we have neglected the direct effect of the oceanic 
drag on the wind field ( F a = 0) as frictional effects 
are included in the internal linear drag term. Also, 
because the effect of the atmospheric heat flux on 
the ocean cannot be properly represented without 
some kind of explicit oceanic mixed-layer dynamics, 
we have assumed Q 0 = 0 and concentrate only on 
the wind stress effect. If we assume that the baro- 
tropic component or the external mode of the wind 
in a two-layer model is given by 



<««> = + «a, 



( 8 ) 



and that the wind at the surface is taken to be that of 
the lower layer, then by (3) and (8) we have 



u s = l M(u a ) - u a ). ^ (9) 

The form of the forcing F 0 in (7) implies that the ef- 
fect of the surface wind stress is to act as a body 
force uniformly distributed over the upper layer of 
depth Dj. Assuming the change in temperature of 
the ocean is much larger in the upper than in the 
lower layer, the first term in the right-hand side of the 
heating function Q a represents a measure of the heat 
flux from the ocean surface. We further assume 
that the perturbation heat flux H n from the ocean is 
given by the form 



h, = Pa c D c„\v,\e u (io) 
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where C„ is the specific heat of air. Using (6) and 
(10) and the hydrostatic relation, the coupling con- 
stant or- will be given approximately by 



2gftp a C 0 |V,| 

p,aD l 



(ID 



where p , is the surface pressure and R the gas con- 
stant. The second term Q c ' in the heating function 
Q a represents the effect of diabatic heating from 
cumulus convection. Thus the form of the forcing 
functions (7) gives a system where the ocean is driven 
by the wind stress and the atmosphere by heat flux 
from the ocean and heating from cumulus convection. 

Next, we determine the projection coefficients 
and 2Pg which arise as a result of the different 
meridional scales of motions in the ocean and in the 
atmosphere. The meridional structure of the uncou- 
pled Kelvin waves in an equatorial p plane has 
the form 



(w fl ,0) ~~ exp(-^y a 2 ) 
(u 0 ,h) - t\p(- l Ay 0 2 ) ’ 



where y a = p l,2 (gH a r il *y, y 0 = p m {gH 0 )-"*y . 

When the ocean forces the atmosphere, the ap- 
propriate projection is 




exp {- l /iy 2 ) exp(-y 2 /2y) dy 




(13) 



where y = ( H 0 /H a ) 1/2 is the ratio of the uncoupled 
phase speeds of the oceanic to the atmospheric 
Kelvin waves. 

Similarly, when the atmosphere forces the ocean, 
the coefficient is 



S = TT~ m j + 

. I 2 V ' 2 

1 1 + y) 



exp(-Vi-yy 2 ) exp (-Viy 2 )dy 



(14) 



In the following we seek the solution of the basic 
equations for some simple cases. Eliminating </> and 
« 0 from (1)— (5) yields 




where we have defined t) = gh, c 0 = ( gH 0 ) 



p,C 0 |Y,| 

2p 0 Di 






a 

0 



1/2 



<T 



and F 0 ' = K(u a ). 

Fourier transform of (12), (13) in space and time 
yields 

AX = Y, (16) 

where 



A = ^ -<a> + + c ° 2 * 2 ikcr \ 

\ -ikKc 0 2 (-ioj + D 0 ) 2 + c 0 2 k 2 J 

x = / w a (co,A:) \ / -kQ c '((o,k) \ 

l r)((D,k) ) < \-ic 0 2 kF 0 '(at,k)J 

Solution to the homogeneous part of (14) is given 
by a vanishing determinant of the matrix A, which 
gives a quartic equation in the frequency w of the 
free modes 



[(-/« + D a f + c a 2 k 2 ][(-ico + D 0 ) 2 

+ c 0 2 k 2 ] = (t) 2 k 2 Co, (17) 
where the coupling frequency o> c is defined by 

< d c ~ (/fCr) 1 ' 2 . 

To illustrate the change of the dispersion relation 
with variation in the strength of the coupling we 
consider two separate cases — with dissipation and 
without dissipation. 



a. Zero dissipation 

By putting D 0 = D a = 0, the homogeneous solu- 
tion to (16) can be obtained analytically,. The fre- 
quency equation (17) yields two solutions for the 
phase speed of the coupled modes: 

( C U 1) 2 = , A(Co 2 + Ca 2 ) ± mc<j 2 _ Co 2)2 

+ 4c 0 2 (o c 2 k~ 2 ] 112 . (18) 

The superscripts I and II are used to denote the plus 
and minus roots which correspond to a fast mode and 
a slow mode, respectively, with eigenfrequencies 
given by o> ul = kc ul . As the phase speeds c 0 , c„ are 
both real and provided a» c kc a the phase speeds of 
the coupled modes c‘ and c 11 must also be real. The 
coupling frequency a> c is nonzero only if k, cr are both 
nonzero. From (18) we see that the coupling favors 
the large scale and approaches zero as the zonal 
wavenumber increases. When oj c = 0, the phase 
speeds reduce to those of the uncoupled modes with 
c l = c a ,c 11 = c 0 . Using the values Pa = 1.25 kg m~ 3 , 
y - 10 _I , p 0 = 10 3 kg m -3 , | V* | = 5 m s _l , C D 
= 1.25 x 10- 3 s-\ D, = 500 m and a = 2 x 10^ 
°C _1 , we. obtain k = 2.16 x 10“ 8 and cr = 1.91 
x 10~ 4 s -1 . This results in a typical value of co c 
=(k<t) 112 = 2 x 10“* s -1 . From the typical values of 
c a ~ 15 m s“ l (Madden and Julian, 1972) andc 0 ~ 2 
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Fig. 2. Transit time for the slow (Mode II) coupled waves to 
cross an ocean basin of width L = 10 000 km for different wave- 
numbers k n = mr/L as a function of the coupling frequency 
u) c (10~® s~ l ). The dashed vertical lines represent the values of the 
critical frequency o> c = k n c a on the abscissa. 



m s” 1 (McWilliams and Gent, 1978), and a horizontal 
scale of the width of the equatorial Pacific, the phase 
speed of the coupled modes are c 1 = 15.04 m s" 1 
and c 11 = 1.6 ms" 1 . The phase speed of the fast 
mode c\ therefore, is not very different from those 
of the uncoupled atmospheric mode, whereas that of 
the slow mode c u is significantly reduced from the 
uncoupled values. For different zonal wavenumbers 
n, Fig. 2 shows the times T n = L/c 11 required for the 
slower of these coupled waves to cross an oceanic 
basin of width L = 10 000 km. It can be seen that 
these transit times are substantially prolonged by 
increasing coupling especially at low wavenumbers. 
The effect becomes increasingly prominent as cu c — » 
oj cr = k n c a [k n = n(n/L )] for which cases the cou- 
pling is unrealistically large and T n tends to infinity. 
The frequency c u cr at which c 11 = 0 has yet another 
interpretation. For fixed k and a the condition oc- 
curs when k = (/ccr) 1/2 /c a wherein the free-traveling 
waves become stationary with respect to the earth. 
In the real system the presence of such modes are 
likely to be important in determining the character 
of the resonant response with respect to stationary 
forcings from topography and differential heating. 
From the values of /c, cr and c a given above, the 
length scale of the stationary mode is 8000 km, which 
is of the same order of scale (~~10 4 km) found in 
some observed low-latitude zonal pressure oscilla- 
tions. However, this scale is modified by the pres- 
ence of friction as we shall see in the next section. 
Further increase of cu c beyond k n c a will render the 
system unstable with respect to that particular mode. 
The detailed stability properties of the solution will 
be discussed in the next section. 



b. Nonzero dissipation 

As we have seen in (i), stationary condition (c 11 
= 0) occurs at oj c = a> cr in the inviscid case and we 
shall see in Section 4 later, when this occurs the 
Mode II amplitude becomes infinitely large. In our 
linear analysis finite amplitude at resonance is pos- 
sible only in the presence of frictional dissipation. 
To see the effect of frictional damping on the growth 
rates and frequencies of the eigenmodes of the sys- 
tem (16), we consider the complete solution to (17) 
with nonzero dissipation terms. The values of the 
dissipation coefficients used are D a = 2.3 x 10~ 6 
s _1 and D 0 = 2 x 10 -7 s _l corresponding to dissipa- 
tion time scales of 5 and 60 days for the atmosphere 
and ocean, respectively. The large value of D a can 
be justified in terms of the effect of cumulus friction 
(e.g., Chang, 1977). However the value of D 0 ap- 
propriate for the baroclinic modes in the ocean is 
less certain. As shown below, within specific ranges, 
the behavior of the solution will have the same quali- 
tative features. By introducing a new variable A 
s /cu, Eq. (17) can be transformed into a quartic 
equation with real coefficients, where the complete 
set of roots can be obtained by numerical methods. 
Thus the real part k r of the eigenvalue k gives the 
growth or decay rate and the imaginary part A,- the 
oscillation frequency of the coupled mode. 

The graphs of k r = A r (cu c ), A* = A,(cu c ) for dif- 
ferent wavenumbers are shown in Figs. 3a and 3b. 
The most interesting results are seen in the Mode II 
response (see Fig. 3b). As in the inviscid case the 
eigenfrequencies A, decrease from an uncoupled 
value at <o c = 0 to zero (stationary condition) as cu c 
approaches a critical frequency cu cr and as long as 
cu c ^ cu cr for a fixed wavenumber n , the damping 
time scale A r is independent of cu c . As cu c increases 
beyond cu cr , the solution first bifurcates into two 
stable non-oscillatory modes (A r < 0). One of these 
modes, with A r increasing, soon becomes unstable 
(A r > 0) as cu e further increases;. Comparison of Figs. 
2 and 3b shows that the stable region (cu < cu cr ) is 
increased as a result of the stabilizing effect of fric- 
tional damping. For a coupling time scale of 4-6 
days (cu c ~ 3 - 2 x 10“ 6 s~ l ) the frequency of the 
largest (n = 1) coupled waves is in the range which 
corresponds to a transit time of 75-90 days across 
a longitudinal span of L = 10 000 km. At this cou- 
pling frequency, the effect on the frequencies of the 
waves with n < 3 becomes negligible. For Mode I 
(Fig. 3a) the eigenvalues indicate. damped oscilla- 
tions (A t - =£ 0 and A r < 0) with damping time scale 
and oscillations frequency virtually independent of 
the strength of the coupling. 

In the stationary limit cu = 0, the frequency equa- 
tion (17) can be solved for the stationary wavelength 
A'.,, yielding 
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of eigenvalues for Mode I. Units are 10 6 s“*. 



Fig. 3b. As in Fig. 3a except for Mode II. 




It can be shown from (19) that k is imaginary, com- 
plex or real depending on the magnitude of the cou- 
pling. Fig. 4 shows the variation of k r = R e(k s ) and 
k t = lm{k 3 ) in Mode II as a function of the coupling 
strength oj c . Three functional regimes are immedi- 
ately obvious. They are [also from Eq. (19)] given 
by the following: 

Regime I: a) c /c a < DJc a - D 0 /c 0 . Here k r = 0 
and/:, > 0. The solution represents pure damping of 
a Kelvin wave toward the east (positive x direction) 
as a result of the relatively overwhelming effect of 
dissipation. 

Regime II: DJc a - D 0 /c 0 < c o c /c a < DJc a + D 0 /c 0 . 
In this range, k r > 0 > 0. The solution consists of 
a wave part with an exponentially damped component 
toward the east. As the coupling increases, the rela- 
tive effect of the damping diminishes reflecting the 
increasing importance of the positive feedback (see 
Section 4) effect of the coupling. 



Regime III: a > c /c a > DJc a + D 0 /c 0 . This regime is 
dominated by the positive feedback effect of the 
coupling with k { = 0 and k t > 0. As the coupling in- 
creases the stationary wave is undamped and its 
magnitude is given asymptotically at large values of 
oj c by the nondissipative limit k r = c o c /c a . 




Fig. 4. Zonal scales of the stationary wave as a function of the 
coupling strength. Abscissa is in units of 10“® s“‘. 
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Fig. 5. Time variation of the Fourier coefficients a(t ) for the 
zonal wind and b{t) for the sea level of the longest wave (« = 1). 
(a) Weak coupling (<*>«. = I0 -6 s" 1 ); (b> normal coupling (o> f = 2 
x IO -6 s" 1 ); and (c) strong coupling (o> r = 4 x I0‘ 6 s -1 ). Ordinate 
is in arbitrary units. 



The zonal scale of the wave and its decay length, 
for typical values of the parameters used previously 
are approximately 15 000 and 8000 km, respectively. 
These length scales are in accordance with some 
observed zonal teleconnections in the equatorial 
oceans and atmospheres (Julian and Chervin, 1978). 

4. Initial value problem 

To illustrate some basic properties of the coupled 
system, we consider a simple initial value problem. 



a{t) = 
bit) = 






b(0)ed - - ia( 0) — y 2 

k 



expC-Zco 1 /) + 
exp(-/V/) 



Because the case with nonzero friction does not 
admit simple analytic solution, we consider the prob- 
lem only for the in viscid case. The case with nonzero 
friction will be dealt with numerically in the next 
section. The eigenvectors corresponding to the eigen- 
values oj 1 and <o n in the inviscid case are 



where 




and Z" = 




S UI = l A(c a 2 - c 0 2 )* 2 



=5= V4[(c« 2 - Co 2 ) 2 * 4 + 4c 0 2 * 2 a) c 2 ] u2 . (20) 

The general solution to the free mode of (16) in 
terms of a truncated Fourier series have the form 




where 

k n = niriL , 

( w*!) = c 'n x Q exp (-i(o'„r) 

+ C" x Z” exp(-/o)||/). (22) 

In (22), C l n and CJJ are complex constants which 
can be determined from initial conditions. To see 
the relative amplitudes of the two modes and their 
temporal variations, we consider the simple case in 
which the coupled system is initially in steady state, 
i.e., 

M0) = a n ( 0) = 0 (23) 

Fig. 5 shows the time variations of the Fourier 
coefficients a(t) y b(t ) of the longest wave ( n = 1) 
for different values of cu c . The most noticeable fea- 
tures are the two-time scale variation in the surface 
wind, a(t) and the Mode II dominated slower varia- 
tion of the sea level b(t). Comparing the three cases 
in Fig. 5 it can be seen that as the coupling is in- 
creased the period of the Mode II response is pro- 
longed and its amplitude becomes larger. 

Next, we consider the structure of the coupled 
waves at the stationary limit. This occurs for strong 
coupling or small wavenumbers. Defining f 
= (o c (kc a )~\ we have from (18) and (20) 

to 1 -» kc a , (o n -» kc 0 ( 1 - f 2 ) 1 ' 2 1 
S'^-cJk 2 ?, S" -* k 2 c a 2 

The solution remains stable provided f < 1. In the 
stationary wave limit (£ —* 1) a>" — > 0, it can be 
shown that (22) can be written in the form 



-a(0)y 



(1 - £ 2 ) ,/2 

K 



b( 0) + ia( 0) 



ikcrb(0) 

k 2 c a 2 

y 



k (1 - $ 2 ) 



1(2 



exp(-/&> u 0 

exp(-/a>"/) 



• (25) 
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Fig. 6. Schematic of the coupling mechanism between the atmospheric and 
oceanic Kelvin waves in the Mode II regime. 



It is clear from (25) that the most important effect of 
the increasing coupling on the ocean is to increase 
the time scale of the purely oceanic mode and to 
enhance the amplitude of the Mode II response by a 
factor of (1 - f 2 ) -1 ' 2 . Resonance conditions occur at 
the Mode II response both in the wind field a(t) 
and the sea level b(t) at £ = 1. Since the fast mode 
(Mode I) in the oceanic response b(t) is of the order 
^0 (y 2 ), they only provide a negligible contribution 
of not more than 1% of the total variation. Physi- 
cally, this is reasonable to expect because the large 
inertia of the ocean compared with that of the at- 
mosphere (y~ ! > 1) makes it unlikely for the fast 
mode to develop there. 

To see the dependence of the coupling on the scale 
of the disturbance in the limit f ► 1 it can be shown 
using the definition of f and (25) that 




!(k 2 c 2 - co c 2 )- 1 , 2 /r l \ 

l i{k 2 c 2 - co 2 )-" 2 I 



exp(-co u r). 



(26) 



which means that for a fixed coupling frequency 
<D c (<kc a ) stronger amplitude is favored for smaller 
values of wavenumber. 

A further result of the coupling is to introduce an 
additional phase factor in the Mode II time variation 
of b(t) given by 



<Pb 



tan -1 



a(0)Ky 

b(0)k(l - £ 2 ) 112 . 



(27) 



On the other hand, the phase of the Mode II re- 
sponse of a(t) is given by 



<Pa 



tan -1 



kcr 



Moxi - ev 2 m 

a(0)yf 2 ' . 



(28) 



It is interesting to note using (27), (28) and the defini- 
tions of f and (d c that the phase difference between 
the dominant waves in the two media is a constant 
given by 

<Pa - <Pb = tt/2, (29) 



i.e., in the .Mode II regime the oceanic thermocline 
response always lags that of the surface wind by 
one-fourth of a wavelength. 

This condition for the phase difference of the Mode 
II waves can be shown to be a very good approxi- 



mation for all ranges of c o c and k which give stable 
solutions. The positive feedback between the oceanic 
and atmosphere Kelvin waves implied by (29) can be 
illustrated as follows. Fig. 6 shows the structure of 
the coupled waves in the Mode II regime. In the 
ocean a high ( H ) will denote a high sea level (baro- 
clinic pressure) or high SST. Since these are Kelvin 
waves the pressure gradient perpendicular to the 
equator are in geostrophic balance. The upper layer 
flow therefore is eastward as a result of which warm 
water is removed from the west of the high and built 
up east of it, causing the wave to move eastward. 
Since the dominant mode has the atmospheric wave 
l A wavelength ahead of its oceanic counterpart it is 
always in a position to assist in the building up of the 
oceanic high which provides increased heat input 
into the atmosphere thus reinforcing its circulation. 
Apparently, the larger amplitude of the reinforced 
oceanic wave makes' it more difficult for the wave 
to be displaced eastward thus increasing the transit 
time of the wave. For a given suitable coupling 
these Kelvin waves in the two media can become 
locked in with the correct phase in a stationary mode. 

5. Numerical experiments 

In this and the following sections we study the 
behavior of the coupled system, i.e., the variations 
of the amplitudes and phase of the coupled waves, 
and see how a model El Nino can evolve given cer- 
tain external forcings. As the. coupled response of 
the model is biased toward ultralong waves and long 
time scales, the spatial domains of the model for 
numerical experimentations are designed along the 
same bias. The large longitudinal extent of the 
response requires that the variables change slowly 
over large east-west distances. As a result they 
should be adequately described by the so-called low- 
resolution models (e.g., McWilliams and Gent, 1978). 
Here we choose the model domains so as to co- 
incide with the bulk regions of the western, central 
and eastern Pacific Ocean (see Fig. 7). Using a 
domain-averaging technique similar to that of Web- 
ster and Lau (1977), the momentum equations in 
(l)-(5) are averaged over each domain and the spa- 
tial derivatives in the height equations are approxi- . 
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Fkj. 7. Schematic representation of model configurations for the 
numerical experiments. 



mated in terms of the domain-averaged quantities. 
The basic equations become 



— + (<f>, - <*>,-,)( A*)- + D a u aJ = 0 

dt 

+ c 0 2 (« a , i+ i - u a ,i)(hx)~ l 

dt 



+ D a <t>i C TTji Qc,i (30) 

+ (Vt ~ tj, + i)(A x)' 1 
dt 

+ Do"o,i + KUa.i- = F 0ti 
— — + C‘o 2 (Wo.l-l-l — ^O.iXA.t) 1 + }i = 0 

dt 

where the overbars denote domain-averaged quan- 
tities, A* is the longitudinal extent of each domain 
and the i index for tj and <b are staggered with respect 
to the domain index denoted by Roman numerals. 
We have taken into account the conditions u a = u 0 
= 0 at the boundary by setting u a%l to = -u a , 2 , u a , e 
= 5, w 0 ,i = -W0.2 and w 0 .2 = -«o,5 on two imagi- 

nary domains outside the boundary. Although the 
coupled equations cannot model explicitly the re- 
flection of the equatorially trapped waves at the 
boundaries, we can still mimic the time scales of 
such process by including the following conditions 
on the phase speed such that 



in accordance with the observation that the oceanic 
Kelvin waves are associated with eastward propa- 
gation, rising and falling sea level in the East and 
West Pacific, respectively, and that the phase speed 
of the gravest Rossby modes is approximately 
one-third of that of the Kelvin waves (Cane and 
Sarachik, 1976). 

In the tropical atmosphere, an important source 
of energy that maintains the east-west overturning 
is in the diabatic heating from latent heat release in 
cumulus convections. Therefore, the interaction pic- 
ture will not be complete without including this 
important effect. Here we assume the following 
simple form of condensation heating, viz ., 



Qr.i — 



0, 



8 f 2? 0 
8^0 



i = 1, 5 



.(32) 



where 5, = (w, Jti+1 - u ati ) is a measure of the upper 
level divergence between adjacent domains and e a 
scaling constant, chosen to be 3 x 10" 5 m s -2 which 
gives a heating rate comparable to that due to adia- 
batic effects. The parameterization in (32) will crudely 
take into account the convective heating in the rising 
branch (positive divergence aloft) of the overturning 
motions. 

Experiments were performed by numerical inte- 
gration of the coupled ordinary differential equations 
(30) with a time step of I day for the atmosphere and 
10 days for the ocean. In the following sections we 
present results for two experiments to study 1) the 
response of the model to sudden relaxation of im- 
posed surface wind stress, 2) sensitivity of the model 
to changes in various parameters. 



6. Results 

In the past, a number of ocean spin down experi- 
ments (e.g., Godfrey, 1975; Hurlbert et al . , 1976) 
have been done using fixed atmospheric boundary 
conditions. The following experiments, however, 
are designed to elucidate the role of the coupling 
mechanism in the light of the development of the 
previous sections. They are not aiming at any 
realistic simulation. Such effort will require more 
sophisticated model structure and physics. For our 
purpose on El Nino will mean an anomalously high 
SST or equivalently a high sea level in the eastern 
Pacific (domain IV). 



a. Wind relaxation 



Co “* 



c K (=2 m s" 1 ) if ^ < 0 
dt 

and > 0 
dt 



Vscu 



otherwise 



In this experiment the coupled model is first spun 
up with constant mean easterly surface wind ( u a ) 
= -6.4 m s*" 1 in the western domains (I and II) and 
{w a } = -4 m s" 1 in the eastern domains (III and 
IV). These values are obtained as averages from the 
25-year surface wind stress of Wyrtki (1975). After 
the steady state is reached the spindown phase for 
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Fig, 8a. Longitude-time section of the sea level during spin- 
down for INT. The time axis is in units of years with year 0 
denoting the start of the wind relaxation. Contours are in 5 cm 
intervals with negative deviations shaded. 

the interactive ocean-atmosphere system (INT) is 
initiated from the steady-state conditions by suddenly 
turning down the mean wind over the western por- 
tions of the Pacific (domains I and II) to -2 m s _1 . 
The same procedure is then repeated with a non- 
interactive atmosphere (NINT) kept fixed at the 
steady-state conditions. At the end of three simulated 
years in the spinup phase a steady state is reached 
with a strong westward equatorial current (—30 cm 
s" 1 ) in the central Pacific and generally increasing 
thermocline depth toward the west. The atmosphere 
is dominated by a single Walker cell with strong 
upper level westerlies in the central Pacific, rising 
motion in the west and sinking motion to the east. 
Figs. 8a and 8b show the longitude-time section of 
the sea level and surface current during the spin- 
down phase. Immediately following the wind relaxa- 
tion in INT (Fig. 8a) a rapid equalization of the 
east-west pressure gradient occurs as a Kelvin wave 
progresses toward the east and eventually causes 
the sea level gradient and the current to reverse 
(Fig. 8b). A seiche-like oscillation follows and be- 
comes strongly damped at the end of the third year 
after spindown. The model El Nino is marked by a 
rapid rise (fall) in sea level in the east (west) fol- 
lowed by a less rapid reversal reflecting the larger 
phase speed of the Kelvin waves compared to those 
of the Rossby waves. Here the triggering of oceanic 
Kelvin waves causes rapid reversal of pressure gra- 
dient and velocity field in about 2-3 months. Hurl- 
bert et al. ( 1976), however, found in a similar oceanic 
spindown experiment that current reversal near the 
equator occurred in only —50 days. The difference 
is most likely due to the interactive role of the at- 

i4- 
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Fig. 8b. As in Fig. 8a except for the surface current. Contour 
interval is 10 cm s' 1 with the eastward flow shaded. 

mospheric Walker-type circulations whereby the 
wave speed is reduced. The times for the coupled 
waves to be effective in the basin are approximately 
3.5 and 10.5 months, respectively, for the Kelvin 
and the Rossby waves giving a periodicity of —14 
months for N the coupled sloshing motion. Observa- 
tions in a number of El Nino's in the past (Wyrtki, 
1975) showed that anomalous conditions were re- 
peated in succession — 12- 15 months apart, and that 
the onset of SST and sea level anomalies, in general, 
occurred faster than its decay over a large extent 
of the equatorial Pacific. Also, of interest is the re- 
sponse of the interactive atmosphere (Fig. 8c). Dur- 
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Fig. 8c. As in Fig. 8a except for the upper level wind. Contour 
interval is 1 ms -1 . Region of easterlies is shaded. 
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Fig. 9a. As in Fig. 8a except for a non-interacting 
atmosphere NINT. 

ing the spin-down the atmosphere is coupled closely 
to the sloshing motions of the underlying oceans — 
an essentially Mode II response superimposed on 
which is a fast (Mode I) response characteristic of 
the atmosphere. In a matter of less than a month 
or so after the maximum thermocline depth in the 
east (El Nino) has been attained, the upper level 
zonal winds are all reversed to easterlies, signaling 
a weakening or even a reversal of the Walker cell. 
The change in wind over domain II is m s' 1 which 
agrees somewhat with the observed averaged 200 mb 
'wind change of —6 m s' 1 over Canton Island during 
periodic anomalous warm waters in the eastern Pa- 
cific (Julian and Chervin, 1978). At the end of the 
integration period, the ocean/atmosphere settles 
down to a steady state with much weakened wind 
and current fields. An important difference between 
the final and the initial steady state in the atmosphere 
is in the extent and the apparent eastward displace- 
ment of the upper level westerlies. This displace- 
ment will imply an eastward shift of the ascending 
branch of the Walker cell whose weakened intensity 
is then consistent with the reduced large-scale east- 
west gradient of SST. The effect of the Walker cir- 
culation to accentuate the initial perturbation in INT 
is evident by its much larger response amplitudes 
compared to NINT (Figs. 9a and 9b). Here the model 
El Nino does not produce the second peak in sea 
level and surface current as in INT. By comparing 
Figs. 8 and 9 it can also be seen that larger spatial 
scales, in addition to longer time-scale oscillations 
are preferred by the interactive mode. In the present 
coupled system positive feedback between the large- 
scale atmospheric and the oceanic waves is made 
possible by the exchange of heat and momentum 



at the air-sea interface. A slackening of surface 
easterlies in the west will cause h to decrease (in- 
crease) in the west (east), decreasing the zonal SST 
contrast which in turn weakens the Walker cell. As 
the SST continues to increase in the central and 
eastern domains, the zonally oriented overturning 
motions become centered over the central and east- 
ern Pacific, thus further diminish the surface east- 
erlies in the west. 

It can be seen in this experiment, the large-scale 
zonal SST gradient is the primary driving force of the 
Walker cell which in turn enhances the oceanic 
response and prolongs its time scales. The key point 
to note here is that in this coupled system an oceanic 
Kelvin wave has to be first triggered before the posi- 
tive air-sea feedback can enhance its amplitude. In 
our model the triggering is achieved by relaxation 
from a period of sustained easterlies in the baro- 
tropic component of the wind. While this suggests 
that the ocean anomaly is caused by the atmosphere, 
it does not contradict the fact that the Walker cell is 
controlled by the large-scale SST gradient. What is 
manifest in this experiment is the strong baroclinic 
response of the coupled system in the presence of 
positive feedback as a result of a change of its baro- 
tropic wind component. The reversal of the zonal 
wind after the maximum sea level in the east also 
agrees with the observational study of Newell and 
Weare (1976) which showed that in some cases at- 
mospheric changes are preceded by SST changes in 
the equatorial region. 

b. Parameter tendencies 

There are eight important parameters used in the 
model viz c a , c 0 , D a , D 0 . C D , k, cr and €. Typical 
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Fig. 9b. As in Fig. 8b except for NINT. 
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values of the first five of these parameters can be 
obtained from the literature and have been listed in 
Section 3. Only the last three are new parameters. 
The largest uncertainties in k and cr according to 
their definitions (Section 3) are in the uncertainty of 
determining the depth of the wind-driven upper layer 
of the ocean. However, their values are restricted to 
within 20% of the values used to give a realistic 
coupling time scale of the order of months. The 
parameter e governing the convective heating in (32) 
is the only free parameter in the model. It is therefore 
desirable to determine the sensitivity of the model 
with respect to this parameter. A series of nine 
separate model runs are made each similar to that 
described in Section 6a. In each of these runs the 
coupled model is spun up from rest with the same 
wind forcing as in Section 6a and then relaxed after 
steady state is reached with the condensation 
heating parameter e = 0, 6 0 , 2e 0 (e 0 = 3 x 10' 5 m s' 2 ) 
under conditions of zero (co c = 0), normal (cu c = 2 
x 10 -6 s' 1 ) and strong (to c = 4 x 10' 6 s' 1 ) coupling. 
The sensitivity of the model with respect to this 
parameter is evaluated in terms of the amplitudes of 
the sea level changes and atmospheric wind changes 
between each generated model El Nino and its 
corresponding initial steady state. Given that the 
coupling mechanisms and the convective heating 
parameterization are valid approximations, the 
result of these experiments will provide further 
insight into the relative importance of these 
processes in the real system. 

Tables la and lb show the amplitude of the changes 
in sea level A h and zonal wind A u a during the Model 
El Nino. For the wind changes, it should be remem- 
bered that A u a is the change only of the baroclinic 
component. The total change is obtained as the av- 
erage of A u a and the change in the barotropic com- ' 
ponent A <w a ), i.e., -4.4 m s' 1 in domains I and II 
and zero in domains III and IV. In all cases the 
distribution shows that the sea level drops in the 
west and rises in the east and that A u a is everywhere 



Table la- Anomalous sea level changes (cm) at different east- 
west locations of the model El Nino for different convective 
heating parameters € and different degrees of air- sea coupling. 
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Table lb. Same as in Table la except for the zonal wind (m s" 1 )* 
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negative, indicating a weakening or reversal of the 
Walker circulation with decreased (increased) east- 
erlies below (aloft). The maximum sea level drop is 
in the west as a result of coupled oceanic Kelvin 
waves leaving the western boundary and the maxi- 
mum wind change is over the western central Pacific. 
This pattern is not changed by a ± 100 % change of e 
from its reference values e 0 . It is clear that the 
amplitudes of the oscillations are enhanced by 
increased coupling and by increased condensation 
heating. The blanks in the tables for the case of 
strong coupling [<o c = 2co ( c 0) ] and strong heating 
(e = 2€ 0 ) correspond to an unstable situation and 
suggest that the model is more sensitive at these 
limits. As shown in Section 6a, the convective heat- 
ing also causes wind fluctuations in the fast time 
scale. Thus the heating contributes both to the 
amplitudes of Modes I and II. The Mode II 
contribution, however, is only possible in the 
presence of coupling. For the case of zero coupling, 
there is no change in the baroclinic component of the 
atmosphere and the amplitude of the sea level 
oscillation is much smaller. The case with normal 
coupling and heating gives a realistic change of 
A u a = -5 m s' 1 in the western central pacific 
(domain II) and A h ~ -40 cm in the extreme west. 
The latter will correspond to a SST change of ap- 
proximately -4.5°C over upper layer ocean depth of 
400 m — a typical amplitude of the anomalous SST 
change in El Nino. 

It is important to note that in this model there are 
two sources contributing to the diabatic heating of 
the atmosphere, viz*, the condensation heating given 
by (32) and the heat input directly from the ocean 
which is proportional to the SST. The foregoing 
analyses show that the former contributes both to 
the fast and slow responses whereas the latter affects 
essentially only the slow time scale because it is 
related to an oceanic variable. However, physically 
the distinction of the two forms of heating is not 
mutually exclusive, as the heat input from ocean to 
atmosphere has to x be realized in some type of 
condensation processes. 
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7, Summary and concluding remarks 

In this paper we have studied the response of a 
simple two-component ocean/atmosphere system 
applicable in the equatorial regions to free and forced 
conditions. The heuristic model employs four basic 
equations, governing the zonal momentum and the 
temperature variation for the atmosphere and the 
ocean, respectively. Within each component, the 
equations are the shallow-water equations with sim- 
plifying assumptions based on the observed proper- 
ties of large-scale equatorially trapped waves in 
relation to short-term climatic time scales. In prac- 
tice the model can be compared to the internal mode 
of a two-level baroclinic system of equations within 
each component. The equivalent depths are chosen 
to correspond to the dominant baroclinic modes in 
the atmosphere and ocean respectively. In the ocean, 
provided the lower layer has large stability, the 
solutions can be regarded as vertical standing modes 
(Philander, 1978). In the atmosphere Chang (1976) 
showed that tropospheric forced wave has no verti- 
cal propagation within the troposphere. In the pres- 
ence of frictional damping a mode exists with large 
vertical wavelength and strong vertical trapping 
(Chang, 1977). Therefore, the present atmospheric 
solution can be identified with this viscous mode. 
The energy exchange between these oceanic and 
atmospheric modes are then effected by the coupling 
in surface wind stress and heat flux. 

Our analyses show that the inclusion of air-sea 
coupling into the linearized shallow- water equations 
can result in two types of dispersion relations for 
the Kelvin waves. The first type (Mode I) has fast 
phase speed and is mostly manifest in the atmos- 
pheric response, being relatively unaffected by air- 
sea coupling. The second type (Mode II) has slower 
phase speed and is the predominant mode in the 
equatorial climatic time-scale. The time-scale and 
the amplitude of the latter mode is substantially in- 
creased by increasing air-sea interaction. Because 
of the widely differing response time of the ocean 
and the atmosphere the effect of the coupling is to 
superimpose on the fast variation of the latter, a 
variation of a longer time scale relevant to oceanic 
processes. Most important is the existence of a 
positive feedback between the atmospheric zonal 
circulation and the large-scale ocean circulation 
in the equatorial plane, leading to the resonant 
excitation of a stationary mode within the coupled 
system. It is plausible that this mode is related to 
the quasi-stationarity of low-latitude teleconnec- 
tions in the real system. 

The results of the numerical experiments using 
the simple coupled model show that an El Nino 
type of oscillation can occur in a baroclinic ocean- 
atmosphere system in the form of trapped equatorial 
waves as a result of a prolonged period of weakening 



in the barotropic component of the wind stress. 
Specifically, the antecedent strengthening and en- 
suing relaxation of surface easterlies in the central 
and western Pacific and the subsequent triggering 
of a coupled Kelvin wave lead to a rapid onset of 
El Nino conditions in the equatorial Pacific. This 
is a simple verification of Wyrtki's (1975) hypothe- 
sis. The large amplitude in the oscillation is a result 
of the coupling through exchanges of heat and mo- 
mentum from the ocean surface. The weakening and 
eastward shift of the rising branch of the atmos- 
pheric Walker cell during El Nino provides a posi- 
tive feedback favoring warm water formation in the 
eastern Pacific. In the experiment, the warm SST 
development is checked by the reflection at the east- 
ern boundary of westward propagating equatorial 
trapped Rossby waves. In practice, the warm water 
likely will be dispersed poleward as well as west- 
ward along the equator by a combination of Rossby 
and Yanai wave modes (Godfrey, 1975; Cane and 
Sarachik, 1976). 

Results also showed that the basic driving force of 
the equatorial Walker cell is the differential diabatic 
heating which is derived mostly from the largescale 
SST contrast. The recent theoretical study of the 
zonally symmetric circulation in a coupled airsea 
system by Lau (1979) and the observational study by 
Comejo-Garrido and Stone (1977) showed that 
indeed there is a high correlation between condensa- 
tion heating and SST. However, the present findings 
differ significantly- from those of Comejo-Garrido 
and Stone (1977, hereafter referred to as CGS) to 
warrant a more careful discussion. CGS indicated 
that while SST may be important in determining the 
initial response of the Walker circulation, it only 
plays a minor role in maintaining the circulation. 
They suggested that the release of latent heat is the 
primary driving force of the Walker cell, but the 
source of the heating is derived not directly from 
evaporation from the ocean, but largely from dy- 
namical moisture convergence. As indicated above, 
evaporation and dynamical convergence are included 
in our model. Both are found to be driving forces for 
the Walker cell. Yet in the time scale of interest, i.e., 
the oceanic time scale, our model shows that the 
circulation is controlled strongly by the heat 
contribution directly from the surface evaporation 
which is proportional to SST. The two apparently 
contradictory views can be resolved by the fact that 
the mechanisms discussed in this paper allow 
transient structure to exist such that the forcing 
functions, i.e., the wind stress and the oceanic heat 
flux are functions of the phenomenon itself and not 
merely boundary conditions. The heat budget study 
of CGS assumed time-averaged conditions and 
therefore is expected to give only steady-state results. 
In the time scale of our study the oceanic time scale 
becomes important and the stationarity assumption 
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is no longer valid. However, as stated by CGS, due 
to the presence of cloud cover, the evaporation over 
regions of high SST can be significantly reduced, 
thereby reducing the effectiveness of the coupling 
and our result may be significantly changed by in- 
cluding this effect. Nevertheless, such inclusion re- 
quires a more detailed radiation calculation in order 
that the effect of solar and infrared radiations can be 
handled properly. In such a case a compatible ocean 
model should also consist of more detailed mixed- 
layer dynamics. The inclusion of these physical 
processes into a coupled model can perhaps shed 
light on some features not predicted by the present 
model, e.g., the observation that heat storage and 
evaporation in the eastern equatorial Pacific did not 
increase during the 1972 El Nino (e.g., Wyrtki et al . , 
1977). Another obvious weakness in this model is the 
absence of an explicit moisture budget thus 
neglecting the important effect of dynamic moisture 
convergence. 

While the results of the analyses described in this 
paper illustrate the causal sequence involved in the 
purely equatorial and zonal processes of the ocean/ 
atmosphere in a very simple manner, because of 
the limited model physics, the ability of the model to 
simulate the real system is seriously deficient. Never- 
theless, the simple model may provide some useful 
backgrounds for data interpretation and for develop- 
ing more sophisticated models of the ocean-atmos- 
phere system. 
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ABSTRACT 

The linearized shallow-water equatorial /9-plane equation was solved for a subset of approximate solutions 
applicable to thermally driven large-scale tropical circulation. In particular, the heat-induced monsoon 
circulations during Southeast Asian northeasterly cold surges are investigated. It was demonstrated that 
the response of the tropical atmosphere to a localized heat source consists of forced Rossby waves propagating 
westward and Kelvin waves eastward along the equator away from the region of forcing. In general, for 
any source/sink distribution, the heat-induced motion can have the characteristics of a Walker-type ( v 
= 0 at the equator) or a Hadley-type (u = 0 at the equator) response or a combination of both, depending 
on the latitudinal location of the forcing. Away from the equator, a forcing corresponding to the sudden 
imposition of mass at the lower layer, or equivalently in our model a rapid cooling of the lower troposphere, 
produces a sudden local surface pressure rise and strong anticyclonic flow to the west of the forcing. Strong 
NE-SW tilt in the axis of the anticyclone is observed and can be understood in terms of the dispersion of 
the various wave modes excited. The low-latitude response is, as expected, dominated by Kelvin and the 
gravest Rossby wave modes. 

Coupled with an equatorial heat source, the sudden cooling of the lower troposphere over a localized area 
in the subtropics gives rise to a northeasterly wind surge and large-scale Walker and Hadley circulations 
reminiscent of periods of strong cold surges over East Asia. Finally, the effect of the presence of a mean 
wind is shown to modify the spatial extent of the equatorial circulation with a mean easterly wind favoring 
the formation of equatorially trapped Walker cells. 



1. Introduction 

It has long been known that a large part of the 
tropical circulation is thermally driven. The observed 
large-scale quasi-stationary zonally asymmetric mo- 
tions derive most of their energy from diabatic heat- 
ing processes within the tropics (Ramage, 1968; 
Webster, 1972). Major condensation heat sources 
within the tropics are found to be over the three 
continents straddling the equator, i.e.. South Amer- 
ica, equatorial Africa and the “maritime continent"’ 
of Borneo and Indonesia (Krishnamurti et al. y 1973; 
Krueger and Winston 1973, 1975). Although these 
heat sources are relatively confined spatially, the 
motions that they drive are of planetary scale. Lon- 
gitudinally localized meridional circulations (local 
Hadley cells) extending to midlatitudes and equa- 
torially trapped east-west oriented circulations 
(Walker cells) covering half of the globe emanate 
from these source regions and constitute a major part 
of the tropical general circulation. Since the zonal 
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asymmetries of these circulations arise as a result of 
land-sea thermal contrasts which determine the heat 
source/sink distribution, these heat-induced motions 
can essentially be regarded as part of the planetary- 
scale monsoon circulation. In particular, during the 
Northern Hemisphere winter, the diabatic heat 
source over the Borneo-Indonesia region is by far 
the strongest of the three tropical heat sources (Ra- 
mage, 1968). The large thermal contrast between 
this equatorial heat source and the immense cold 
source over northern China and Siberia provides the 
major driving force for the East Asia winter mon- 
soon, which is one of the most energetic wintertime 
circulation systems in the tropics. The understanding 
of the dynamics of this circulation is crucial to fore- 
casting the weather in this part of the world. 

Although observational work on tropical large- 
scale circulations is abundant (e.g., Krishnamurti et 
al. y 1973; Krueger and Winston, 1973, 1975; Mu- 
rakami and Unninayar, 1977; Murakami, 1980a,b; 
and many others), theoretical studies aiming at dy- 
namical interpretation of these circulations are rel- 
atively few. The pioneer work of Webster (1972) 
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suggested that the Walker circulation may be re- 
garded as rotationally trapped Kelvin waves in re- 
sponse to diabatic heating. However, Webster did 
not provide an explanation for the observed east- west 
asymmetry in his results. Further, the vertical struc- 
ture of his Kelvin wave did not agree well with ob- 
servation. Chang (1977) showed that in the presence 
of frictional damping and Newtonian cooling, the 
forced Kelvin wave indeed has the correct phase 
speed and vertical wavelength. In a recent study Gill 
(1980) addressed the zonal asymmetry problem by 
pointing out that the steady state response of the 
tropical atmosphere to an isolated heating function 
concentrated at the equator consists of planetary 
Rossby waves to the west and Kelvin waves to the 
east of the forcing. All these studies point to the fact 
that at least to a first approximation, the tropical 
large-scale circulation, including the monsoons, can 
be regarded as being driven by local heat sinks and 
sources, provided they are related realistically to the 
actual land-sea distribution. However, for more de- 
tailed and accurate descriptions the effects of orog- 
raphy and lateral coupling with midlatitudes will 
have to be included. 

In this paper we seek simple solutions appropriate 
to thermally-forced long waves in a divergent bar- 
otropic equation in an equatorial /3-plane. Here, we 
stress the planetary-scale aspect of the heat-induced 
motions with a particular solution relevant to thermal 
forcings during the winter monsoon period. As the 
model is highly idealized, we do not attempt to make 
detailed simulation of every nuance of the monsoon 
circulation, but rather we aim at the physical inter- 
pretation of the large-scale dynamical processes in- 
volved. In Sections 2 and 3 the model is developed 
and general solutions applicable to all heating func- 
tions are obtained. Section 4 focuses on the char- 
acteristics of the response in the Hadley and the 
Walker regimes. In Sections 5 and 6, we start with 
a brief discussion of the East Asian monsoon and its 
relevance to thermally-forced circulations. Solutions 
using specific forcing, functions are then obtained and 
results analyzed. In Section 7 we examine the crucial 
effect of a mean zonal wind on the structure and 
intensity of the Walker circulation. 



2. The model 



The basic equations are the shallow-water equa- 
tions in an equatorial /3-plane: 



du rr du „ d<b 

- + U--i3yv + - + 0, 

dv ;,dv d<t> 

- + 0, 

d<i) d<t> ( du dv\ 

i +u T x +gh [y x + ^) +i '^ Q ' 



(1) 

( 2 ) 
(3) 



where u, v are the zonal and meridional velocity lin- 
earized from a basic state with constant zonal wind 
[/; <f> is the geopotential, €, e f the Rayleigh friction 
and Newtonian cooling coefficient, respectively; /3 
= 2 Q/fl, where Q is the angular frequency and a the 
radius of the earth; Q the forcing (mass sink or 
source) term, and h an equivalent depth for the trop- 
ical atmosphere. Although (l)-(3) are the so-called 
divergent barotropic equations, their validity as 
representing a particular baroclinic mode corre- 
sponding to a specific structure in a continuous strat- 
ified atmosphere or as a two-layer atmospheric 
system has been demonstrated by Matsuno (1966) 
and Gill (1980). Here we prefer the two-layer in- 
terpretation, whereby heating means increasing the 
mass of high potential temperature fluids or equiv- 
alently transferring mass from a lower to an upper 
layer, and the reverse situation occurs for cooling. 
Because of the limited degree of freedom in the shal- 
low-water system, the implied upper and lower flows 
have to be exactly opposite in order to maintain con- 
tinuity. We should point out that, in the subsequent 
discussion, when a mass source is referred to in the 
surface layer, cooling of the troposphere is implied, 
and when it is referred to in the upper layer, heating 
is implied. 

The eigenfrequencies co of the free modes of the 
system ( 1 )— (3) are given by the dispersion relation 
for 

a; = (w r — Uk) + Zoj, , (4) 

as 

(oj, - Ukf 
pfgh /3 

+ T^7m' [rh ~ (1 " + ' ) (5) 

and 

oj, = €, 

where k is a wavenumber. 

The dispersion relation (5) for the Doppler-shifted 
frequency (w r — Uk) depicted in Fig. 1 shows that 
for large-scale waves (small k) only the Kelvin wave 
and the Rossby waves have small frequencies. For 
an equivalent depth of 40 m the fundamental phase 
speed of c = tfgh is 20 m s'" 1 , close to that observed 
for Kelvin wave-like oscillations in the tropical tro- 
posphere (Chang, 1977). At length scale greater than 
4000 km, the gravity waves and the mixed Rossby- 
gravity wave all have frequencies corresponding to 
a time scale of less than V 2 day, which is too short 
to be of significance for planetary-scale motions of 
the monsoonal time scale. Therefore if we emphasize 
the slow and large-scale waves (stippled region in 
Fig. 1), the dispersion relation (5) for these waves, 
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Fig. 1. Dispersion relation for the real part of the eigenfrequency 
and the wavelengths in an equatorial /3-plane. The region within 
which the longwave approximation remains valid is indicated by 
stippling. 



to a very good approximation becomes 

Ygh = (In + 1 ). ( 6 ) 

(j) r — UfC 

The “long wave” approximation expressed by (6) 
had been shown (Lighthill, 1969) to be equivalent 
to assuming geostrophic balance in the zonal direc- 
tion, which is a fairly accurate assumption for large- 
scale Rossby- and Kelvin-type waves. We shall show 
in the following that such balance does lead to the 
dispersion relation (6). Following a procedure similar 
to Gill (1980), we non-dimensionalize (l)-(3) in 
terms of the basic equatorial length, time and veloc- 
ity scales 

L = (2 T = Wr xl \gh)-' ,A 



and C = Ygh, 

and define two new non-dimensional quantities, 

i = <f> + u, x = <P ~ (7) 

For large-scale waves and small friction parameter, 
the non-dimensional equivalence of (l)-(3) is the 
following set 

d\b (hi/ dt) 

— + ei + (U + 1) — + x hyv = Q, (8) 

dt dx dy 

dy dy dv ' 

-£ + <X + (U- l)-+ + - + V 2 yv = Q, (9) 
dt dx dy 

j- + Kyi + ~ — ttyx = 0 . ( 10 ) 

dy dy 

where (10) represents the zonal geostrophic condi- 
tion. The non-dimensional frictional parameters for 
wind and geopotential are assumed to be the same 
in (8) and (9). 



To solve (8)— ( 1 0), consider the following expan- 
sions: 

(i, V, Q) = 2 (.in, Xn, V n , Qn)D n (y), ( 11 ) 

n 



where D n (y) are parabolic functions of order n. For 
small values of c, these functions are approximate 
eigenfunctions to the system (8)-(10) and have the 
following properties 

D„(y) = (-1)" exp( x Ay 2 ) ^ exp(-'/z/), (12) 



f D n (y)D m (y)dy = Yl*n\b mn (13) 



and 



dDn 

dy 



+ ViyD„ = /!£>„_, , 



- ViyD n = -D a+l . 
dy 



(14) 



(15) 



Using (11), (14) and (15), after solving for \p„, 
v„ and Xn (8)— ( 10) can be written as 

+ (U + 1) + tio = Qo . (16) 



di, 



n+i 



dt 



(r, 1 . 

( u ■ 5TTT) — + 



_ nQ n +\ + Qn~\ 
(2n + 1) 

*i-0 

2{n + 1)3*^! [( n + \)Q n + x - Q„-,] 



, (17) 



V n = 



2n + 1 dx 
Vo = “01 
and 



2 n + 1 



(18) 



Xn-I = in + 1 )*„+!, for n> 1. (19) 



It is readily seen that (16) represents a Kelvin 
wave moving eastward with unit phase speed and 
(17) represents Rossby waves moving westward with 
phase speed ( 2n + l)' 1 both relative to the mean 
wind U. Notice that in this approximation the Rossby 
wave speeds are independent of the east- west wave- 
number as predicted by (6) and that the high-fre- 
quency gravity modes are completely filtered making 
possible the separation of the large-scale Kelvin and 
Rossby modes in a very simple manner. Given the 
space- and time-dependence of the forcing Q , the 
functions 4 \ y and v can be determined from (16)- 
(19) and the velocity and geopotential fields can then 
be obtained from (7). 



3. The initial value problem 

In this section we consider the general solution to 
a simple switch-on forcing. The solution for arbitrary 
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time-dependence can then be obtained by convolu- 
tion. While the choice of the switch-on time depen- 
dence may not be entirely consistent with our as- 
sumption of slowly varying modes, because of the 
high frequencies contained in such function, we note 
that the “long wave” approximation we used filtered 
all these fast modes or transients in our problem. The 
solutions obtained are only a subset of the total so- 
lution. Nevertheless, this is an important subset, be- 
cause in reality we expect that the largest uncertainty 
in the response of the tropical atmosphere to heat- 
induced forcing is due to the large variability of these 
fast, gravity-like modes arising from even slight var- 
iations in the temporal and spatial distribution of the 
forcing. The degree of sophistication in our present 
model does not warrant a detailed description of 
the these rather complex variations. However, all 
these responses share common characteristics which 
are most likely carried by the slowly-varying modes. 
If we look at the average response over a large num- 
ber of cases with similar but slightly different forc- 
ings, as in composite methods used in observational 
studies, these fast modes would be expected to cancel 
out as “background noise.” Therefore, for the pur- 
pose of this work, the use of the simple time-depen- 
dence in the forcing is not in conflict with our sim- 
plifying approximation used for the solutions. 
Consider a forcing of the form 

Q„(x, t) = q n exp (~a 2 x 2 )H(t), (20) 



where H(t) is the Heaviside step-function given by 

for 
for 



, f 1 for / > 0) 

m o-L , , J, (2D 

10 for t < 0J 



and a" 1 is the east-west scale of the heating function. 
Eqs. (16)-(18) are then easily solved by transform 
methods. The solutions are given by, for the Kelvin 
wave 

u*.t)- * «p{- «7TT)} 

X {(?o(x) - Go[x - (C/ + 1)/]} (22) 



and 



GJLx) = 



' te,p {s(F^)F 



2(2n+ l)a{U-c„) " 

X Erfc {.(« - J • f ” “ » '• < 25 > 



where 



(r 0 = sgn(C/ +1), <r„ = sgn(C/ - c„) 

and 

1 

Cn “ ~ i i ’ 

+ 1 

and Erfc(z) is the complimentary error function de- 
fined by 

Erfc(z) = -p f exp(-/ 2 )^r. (26) 

Vtt Jz 

The cases of (U + 1) = 0 and ( U — c n ) - 0 represent 
stationary resonant situations and require special 
treatment. They are considered in a later section. It 
is clear that for non-resonance the solutions (22) and 
(23) have a stationary component and a transient 
component which move away from the source region 
with non-dimensional phase speed 1 and c n to the 
east and west, respectively. The magnitude of the 
damping'depends on the friction coefficient and the 
propagation speed. For the case of no mean wind 
( U = 0), the ratio of decay distances to the east and 
west for the strongly trapped or the lowest-order 
mode is 3 to 1 , implying strong zonal asymmetry in 
the east-west direction. 

To solve the complete problem for a specific hor- 
izontal heating distribution, the next step is to de- 
termine the expansion coefficients q n in (20). These 
coefficients are most simply obtained for the case 
when the forcing distribution is also Gaussian in the 
^-direction. Consider a general forcing function of 
the form 



and for the Rossby waves 

t) = (nq n+l + <?„_,) expj- ) } 

X {<?„(*) - G„[x -{U- c„)f]}, (23) 

where 



Vtt 



exp 



Go(x) = 



i l 2 

2 a(U + 1 )] 



2 a(U+ 1) 



Q(x, y , t) 

= exp[— (>> - >v) 2 A 2 ] exp (~a 2 x 2 )H(t), (27) 

where y e is the center and A -1 is a measure of the 
width of the Gaussian function in y. The coefficients 
of the heating function in (20) are given by 

1 f +0 ° 

q " = Wlff J- exp ^ - ^ ~ T'c) 2 * 2 ]^ 004V. (28) 

which has a recurrence relationship of the form 



X (X 0 



Erfc] (T 0 - 









(«+!)(!+ 4A 2 )q„ 



X*- 



= 4A 2 ^„ + (1 - 4A 2 )^, , 



(29) 
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with 

<70 = 2 1/2 (1 4- 4X 2 ) exp[— y c 2 X 2 /(4X 2 + 1)] 



and 

< 7 , = 2 1/2 ( 1 + 4X 2 )“ 3/ V c 4X 2 exp[-^ 2 X 2 /(4X 2 + 1)]. 

Given a general forcing function of the form (27), 
or/ a combination of several of these functions with 
different geometric shape factors y c > X, a, the initial 
value problem is completely solved using (22)-(25), 
(28)-(29). In the latter case, the solutions for mul- 
tiple mass sources and sinks can be obtained by linear 
superposition. 



4. Special solutions 
a. Symmetric heating 

The properties of the solutions (22)-(25) are best 
examined by considering, a heating function of the 
following form 

Q{x,y, t) = expt-'/K* 2 + (30) 

which corresponds to a circular heat source concen- 
trated at the equator with radius equal to one equa- 
torial Rossby radius of deformation (y c = 0, X = a 
* Vi). From (29)-(30), we see that all q n 's except 
<7o( - 1 ) are zero. The only non-zero components from 
(22)-(25) and (18) are \p Qy xo» ^ 2 * The geopoten- 
tial and velocity fields are given by 




Kelvin mode 




V 2 (t 2 D 2 (y) + xoDoiy)) 



V 2 (\p 2 D 2 (y) - XoDoiy)) 



V\D x iy) 

Rossby mode (n = 1) 



(31) 



Fig. 2 shows the upper level response to the heating 
at the end of two, five, and nine simulated days. The 
solution approaches the steady state at the end of 
nine days. It is clear from Fig. 2 that the imposed 
heating generates Kelvin and Rossby (n = 1 ) waves 
which emanate to the east and west with different 
phase speeds and markedly different horizontal 
structures. Most noticeable is the strong zonal out- 



flow trapped about the equator with u P v. The u 
and <f> fields are symmetric and the v field antisym- 
metric with respect to the equator. Close to the equa- 
tor the flow is dominated by a Walker-type circu- 
lation. The amplitude of the Hadley-type component 
is relatively small in this case. A different perspective 
of the response of the system can be obtained from 
an examination of the vorticity equation appropriate 
to the system (l)-(3) for U - 0 




Py<t> \ 

gh) 



+ u/3 = S 



m 

gh ’ 



(32) 



where f* is the relative vorticity. It is readily seen that 
for a symmetric heating function with positive Q at 
the upper level (heating), a vorticity sink/source pair 
is produced on opposite sides of the equator. For the 
function (30), these sink and source maxima, ob- 
tained by differentiating the function S with respect 
to y , are at a distance of il times a Rossby radius 
of deformation away from the equator. At the equa- 
tor, however, vorticity generation is zero and the flow 
becomes strongly divergent. These features are rep- 
resented very well in Fig. 2. 



b . Asymmetric heating 

Here we consider the solution for the case in which 
the heating function is of the form 

Q(x,y, t) = jexp[-!40' 2 + x 2 )]/7(/). (33) 

In (33), the ^-dependence is just that of D x (y) from 
(12). By the orthogonality relation (13), the only 
non-zero coefficient is ^i( = 1 ). Using (22), (23), and 
(17), the pressure and velocity fields can be shown 
to have the following form 

f <t>\ f + xi^OO] 

« = VMiDAy) - XxD x {y)] 

\vj \ v 2 D 2 (y) 

Rossby mode ( n = 2) 

( 

X I xpiDi(y) I . 

\ tfoAjOO J 

mixed Rossby-gravity mode (n = 0) 

As VoD 0 = — expf^O' 2 + x 2 )] from (18), the n = 0 
mixed mode in this case represents a purely local 
response. The pressure and zonal wind fluctuation 
of the mixed mode are in fact zero because \p x ^ 0 
by (17). 

Fig. 3 shows the evolution of model response at 
the upper level for days 2, 5, and 9 respectively. Here 
no Kelvin waves are excited and the response consists 
of a local component plus a slowly westward-prop- 
agating component with phase speed c 2 = 0.2 (~4 
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Fig. 2. Upper-level wind and geopotential response to a heating function of the form 
Q - expf-Cx 2 + y 2 )/ 4] at days 2, 5, 9 with the friction coefficient a - 0.5. 



m s” 1 ). The presence of strong cross-equatorial flow 
with anticyclonic centers on the opposite side of the 
equator is eviderft of the antisymmetric nature (with 
respect to u and 4>) of the response. Such features 
are characteristics of strong Hadley-type circula- 



i 

tions. In this case from (32), two upper-level vorticity 
sinks are produced at a distance of twice the Rossby 
radius of deformation on opposite sides of the equa- 
tor, and as the zonal pressure gradient vanishes at 
the equator the zonal flow is zero there. 
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Fig. 3. Same as in Fig. 2 except for a heating function of the form Q - y exp^/H* 2 + y 2 )]H(t). 



Summing up the two cases, it is clear that a sym- 
metric forcing concentrated at the equator will pro- 
duce a symmetric (Walker-type) response. (Kelvin 
waves are excited only if q 0 ^ 0.) Likewise, an an- 
tisymmetric forcing will produce only an antisym- 
metric (Hadley-type) response. For the Walker re- 



gime, symmetric mode Rossby and Kelvin waves are 
responsible for propagating disturbances to the west 
and east of the source region, respectively. Con- 
versely, in the Hadley regime, only the antisym- 
metric Rossby modes are important and an implied 
north-south overturning to the west is set up con- 
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